home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / zsyrk.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  9.1 KB  |  297 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE ZSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA,
  5.      $                   BETA, C, LDC )
  6. *     .. Scalar Arguments ..
  7.       CHARACTER*1        UPLO, TRANS
  8.       INTEGER            N, K, LDA, LDC
  9.       COMPLEX*16         ALPHA, BETA
  10. *     .. Array Arguments ..
  11.       COMPLEX*16         A( LDA, * ), C( LDC, * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  ZSYRK  performs one of the symmetric rank k operations
  18. *
  19. *     C := alpha*A*A' + beta*C,
  20. *
  21. *  or
  22. *
  23. *     C := alpha*A'*A + beta*C,
  24. *
  25. *  where  alpha and beta  are scalars,  C is an  n by n symmetric matrix
  26. *  and  A  is an  n by k  matrix in the first case and a  k by n  matrix
  27. *  in the second case.
  28. *
  29. *  Parameters
  30. *  ==========
  31. *
  32. *  UPLO   - CHARACTER*1.
  33. *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  34. *           triangular  part  of the  array  C  is to be  referenced  as
  35. *           follows:
  36. *
  37. *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
  38. *                                  is to be referenced.
  39. *
  40. *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
  41. *                                  is to be referenced.
  42. *
  43. *           Unchanged on exit.
  44. *
  45. *  TRANS  - CHARACTER*1.
  46. *           On entry,  TRANS  specifies the operation to be performed as
  47. *           follows:
  48. *
  49. *              TRANS = 'N' or 'n'   C := alpha*A*A' + beta*C.
  50. *
  51. *              TRANS = 'T' or 't'   C := alpha*A'*A + beta*C.
  52. *
  53. *           Unchanged on exit.
  54. *
  55. *  N      - INTEGER.
  56. *           On entry,  N specifies the order of the matrix C.  N must be
  57. *           at least zero.
  58. *           Unchanged on exit.
  59. *
  60. *  K      - INTEGER.
  61. *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
  62. *           of  columns   of  the   matrix   A,   and  on   entry   with
  63. *           TRANS = 'T' or 't',  K  specifies  the number of rows of the
  64. *           matrix A.  K must be at least zero.
  65. *           Unchanged on exit.
  66. *
  67. *  ALPHA  - COMPLEX*16      .
  68. *           On entry, ALPHA specifies the scalar alpha.
  69. *           Unchanged on exit.
  70. *
  71. *  A      - COMPLEX*16       array of DIMENSION ( LDA, ka ), where ka is
  72. *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  73. *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  74. *           part of the array  A  must contain the matrix  A,  otherwise
  75. *           the leading  k by n  part of the array  A  must contain  the
  76. *           matrix A.
  77. *           Unchanged on exit.
  78. *
  79. *  LDA    - INTEGER.
  80. *           On entry, LDA specifies the first dimension of A as declared
  81. *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  82. *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
  83. *           be at least  max( 1, k ).
  84. *           Unchanged on exit.
  85. *
  86. *  BETA   - COMPLEX*16      .
  87. *           On entry, BETA specifies the scalar beta.
  88. *           Unchanged on exit.
  89. *
  90. *  C      - COMPLEX*16       array of DIMENSION ( LDC, n ).
  91. *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
  92. *           upper triangular part of the array C must contain the upper
  93. *           triangular part  of the  symmetric matrix  and the strictly
  94. *           lower triangular part of C is not referenced.  On exit, the
  95. *           upper triangular part of the array  C is overwritten by the
  96. *           upper triangular part of the updated matrix.
  97. *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
  98. *           lower triangular part of the array C must contain the lower
  99. *           triangular part  of the  symmetric matrix  and the strictly
  100. *           upper triangular part of C is not referenced.  On exit, the
  101. *           lower triangular part of the array  C is overwritten by the
  102. *           lower triangular part of the updated matrix.
  103. *
  104. *  LDC    - INTEGER.
  105. *           On entry, LDC specifies the first dimension of C as declared
  106. *           in  the  calling  (sub)  program.   LDC  must  be  at  least
  107. *           max( 1, n ).
  108. *           Unchanged on exit.
  109. *
  110. *
  111. *  Level 3 Blas routine.
  112. *
  113. *  -- Written on 8-February-1989.
  114. *     Jack Dongarra, Argonne National Laboratory.
  115. *     Iain Duff, AERE Harwell.
  116. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  117. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  118. *
  119. *
  120. *     .. External Functions ..
  121.       LOGICAL            LSAME
  122.       EXTERNAL           LSAME
  123. *     .. External Subroutines ..
  124.       EXTERNAL           XERBLA
  125. *     .. Intrinsic Functions ..
  126.       INTRINSIC          MAX
  127. *     .. Local Scalars ..
  128.       LOGICAL            UPPER
  129.       INTEGER            I, INFO, J, L, NROWA
  130.       COMPLEX*16         TEMP
  131. *     .. Parameters ..
  132.       COMPLEX*16         ONE
  133.       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
  134.       COMPLEX*16         ZERO
  135.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  136. *     ..
  137. *     .. Executable Statements ..
  138. *
  139. *     Test the input parameters.
  140. *
  141.       IF( LSAME( TRANS, 'N' ) )THEN
  142.          NROWA = N
  143.       ELSE
  144.          NROWA = K
  145.       END IF
  146.       UPPER = LSAME( UPLO, 'U' )
  147. *
  148.       INFO = 0
  149.       IF(      ( .NOT.UPPER               ).AND.
  150.      $         ( .NOT.LSAME( UPLO , 'L' ) )      )THEN
  151.          INFO = 1
  152.       ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
  153.      $         ( .NOT.LSAME( TRANS, 'T' ) )      )THEN
  154.          INFO = 2
  155.       ELSE IF( N  .LT.0               )THEN
  156.          INFO = 3
  157.       ELSE IF( K  .LT.0               )THEN
  158.          INFO = 4
  159.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  160.          INFO = 7
  161.       ELSE IF( LDC.LT.MAX( 1, N     ) )THEN
  162.          INFO = 10
  163.       END IF
  164.       IF( INFO.NE.0 )THEN
  165.          CALL XERBLA( 'ZSYRK ', INFO )
  166.          RETURN
  167.       END IF
  168. *
  169. *     Quick return if possible.
  170. *
  171.       IF( ( N.EQ.0 ).OR.
  172.      $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
  173.      $   RETURN
  174. *
  175. *     And when  alpha.eq.zero.
  176. *
  177.       IF( ALPHA.EQ.ZERO )THEN
  178.          IF( UPPER )THEN
  179.             IF( BETA.EQ.ZERO )THEN
  180.                DO 20, J = 1, N
  181.                   DO 10, I = 1, J
  182.                      C( I, J ) = ZERO
  183.    10             CONTINUE
  184.    20          CONTINUE
  185.             ELSE
  186.                DO 40, J = 1, N
  187.                   DO 30, I = 1, J
  188.                      C( I, J ) = BETA*C( I, J )
  189.    30             CONTINUE
  190.    40          CONTINUE
  191.             END IF
  192.          ELSE
  193.             IF( BETA.EQ.ZERO )THEN
  194.                DO 60, J = 1, N
  195.                   DO 50, I = J, N
  196.                      C( I, J ) = ZERO
  197.    50             CONTINUE
  198.    60          CONTINUE
  199.             ELSE
  200.                DO 80, J = 1, N
  201.                   DO 70, I = J, N
  202.                      C( I, J ) = BETA*C( I, J )
  203.    70             CONTINUE
  204.    80          CONTINUE
  205.             END IF
  206.          END IF
  207.          RETURN
  208.       END IF
  209. *
  210. *     Start the operations.
  211. *
  212.       IF( LSAME( TRANS, 'N' ) )THEN
  213. *
  214. *        Form  C := alpha*A*A' + beta*C.
  215. *
  216.          IF( UPPER )THEN
  217.             DO 130, J = 1, N
  218.                IF( BETA.EQ.ZERO )THEN
  219.                   DO 90, I = 1, J
  220.                      C( I, J ) = ZERO
  221.    90             CONTINUE
  222.                ELSE IF( BETA.NE.ONE )THEN
  223.                   DO 100, I = 1, J
  224.                      C( I, J ) = BETA*C( I, J )
  225.   100             CONTINUE
  226.                END IF
  227.                DO 120, L = 1, K
  228.                   IF( A( J, L ).NE.ZERO )THEN
  229.                      TEMP = ALPHA*A( J, L )
  230.                      DO 110, I = 1, J
  231.                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
  232.   110                CONTINUE
  233.                   END IF
  234.   120          CONTINUE
  235.   130       CONTINUE
  236.          ELSE
  237.             DO 180, J = 1, N
  238.                IF( BETA.EQ.ZERO )THEN
  239.                   DO 140, I = J, N
  240.                      C( I, J ) = ZERO
  241.   140             CONTINUE
  242.                ELSE IF( BETA.NE.ONE )THEN
  243.                   DO 150, I = J, N
  244.                      C( I, J ) = BETA*C( I, J )
  245.   150             CONTINUE
  246.                END IF
  247.                DO 170, L = 1, K
  248.                   IF( A( J, L ).NE.ZERO )THEN
  249.                      TEMP      = ALPHA*A( J, L )
  250.                      DO 160, I = J, N
  251.                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
  252.   160                CONTINUE
  253.                   END IF
  254.   170          CONTINUE
  255.   180       CONTINUE
  256.          END IF
  257.       ELSE
  258. *
  259. *        Form  C := alpha*A'*A + beta*C.
  260. *
  261.          IF( UPPER )THEN
  262.             DO 210, J = 1, N
  263.                DO 200, I = 1, J
  264.                   TEMP = ZERO
  265.                   DO 190, L = 1, K
  266.                      TEMP = TEMP + A( L, I )*A( L, J )
  267.   190             CONTINUE
  268.                   IF( BETA.EQ.ZERO )THEN
  269.                      C( I, J ) = ALPHA*TEMP
  270.                   ELSE
  271.                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  272.                   END IF
  273.   200          CONTINUE
  274.   210       CONTINUE
  275.          ELSE
  276.             DO 240, J = 1, N
  277.                DO 230, I = J, N
  278.                   TEMP = ZERO
  279.                   DO 220, L = 1, K
  280.                      TEMP = TEMP + A( L, I )*A( L, J )
  281.   220             CONTINUE
  282.                   IF( BETA.EQ.ZERO )THEN
  283.                      C( I, J ) = ALPHA*TEMP
  284.                   ELSE
  285.                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  286.                   END IF
  287.   230          CONTINUE
  288.   240       CONTINUE
  289.          END IF
  290.       END IF
  291. *
  292.       RETURN
  293. *
  294. *     End of ZSYRK .
  295. *
  296.       END
  297.